




library(mvtnorm)
library(arm)
library(BRugs)
library(R2OpenBUGS)
library(coda)
library(car)
library(foreign)
library(vioplot)



setwd("~/../Dropbox/iREDS/Data/Analysis/scales/")



# Read in OpenBUGS results and create figure




##  THESE ARE SUBGROUP ONLY
setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S1/wcovariates")
coda.data<-read.openbugs(stem="")
results.S1<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S2/wcovariates")
coda.data<-read.openbugs(stem="")
results.S2<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S3/wcovariates")
coda.data<-read.openbugs(stem="")
results.S3<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S5/wcovariates")
coda.data<-read.openbugs(stem="")
results.S5<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S6/wcovariates")
coda.data<-read.openbugs(stem="")
results.S6<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))

setwd("~/../Dropbox/iREDS/Data/Analysis/scales/S8/wcovariates")
coda.data<-read.openbugs(stem="")
results.S8<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))




setwd("~/../Dropbox/iREDS/Data/Analysis/scales/")


# extract treatment effect estimates and standardize for Cohen's D
std.deltatheta<-c(0.86, 1.17, 1.33, 1.06, 1.16, 1.17)


g.wmpi<- as.data.frame(cbind((results.S1$S1.alpha1)/std.deltatheta[1], (results.S2$S2.alpha1)/std.deltatheta[2], 
                            (results.S3$S3.alpha1)/std.deltatheta[3], (results.S5$S5.alpha1)/std.deltatheta[4], 
                            (results.S6$S6.alpha1)/std.deltatheta[5], (results.S8$S8.alpha1)/std.deltatheta[6]))
names(g.wmpi)<-c("S1.alpha", "S2.alpha", "S3.alpha", "S5.alpha", "S6.alpha", "S8.alpha")



g.urm<- as.data.frame(cbind((results.S1$g.4+results.S1$S1.alpha1)/std.deltatheta[1], (results.S2$g.4+results.S2$S2.alpha1)/std.deltatheta[2], 
                            (results.S3$g.4+results.S3$S3.alpha1)/std.deltatheta[3], (results.S5$g.4+results.S5$S5.alpha1)/std.deltatheta[4], 
                            (results.S6$g.4+results.S6$S6.alpha1)/std.deltatheta[5], (results.S8$g.4+results.S8$S8.alpha1)/std.deltatheta[6]))
names(g.urm)<-c("S1.alpha", "S2.alpha", "S3.alpha", "S5.alpha", "S6.alpha", "S8.alpha")


g.notmale<- as.data.frame(cbind(((-1)*results.S1$g.5+results.S1$S1.alpha1)/std.deltatheta[1], ((-1)*results.S2$g.5+results.S2$S2.alpha1)/std.deltatheta[2], 
                                ((-1)*results.S3$g.5+results.S3$S3.alpha1)/std.deltatheta[3], ((-1)*results.S5$g.5+results.S5$S5.alpha1)/std.deltatheta[4], 
                                ((-1)*results.S6$g.5+results.S6$S6.alpha1)/std.deltatheta[5], ((-1)*results.S8$g.5+results.S8$S8.alpha1)/std.deltatheta[6]))
names(g.notmale)<-c("S1.alpha", "S2.alpha", "S3.alpha", "S5.alpha", "S6.alpha", "S8.alpha")

g.notpi<- as.data.frame(cbind(((-1)*results.S1$g.6+results.S1$S1.alpha1)/std.deltatheta[1], ((-1)*results.S2$g.6+results.S2$S2.alpha1)/std.deltatheta[2], 
                              ((-1)*results.S3$g.6+results.S3$S3.alpha1)/std.deltatheta[3], ((-1)*results.S5$g.6+results.S5$S5.alpha1)/std.deltatheta[4], 
                              ((-1)*results.S6$g.6+results.S6$S6.alpha1)/std.deltatheta[5], ((-1)*results.S8$g.6+results.S8$S8.alpha1)/std.deltatheta[6]))
names(g.notpi)<-c("S1.alpha", "S2.alpha", "S3.alpha", "S5.alpha", "S6.alpha", "S8.alpha")



## Subgroup analysis 

apply(g.urm, 2, mean)
apply(g.urm, 2, sd)

apply(g.notmale, 2, mean)
apply(g.notmale, 2, sd)

apply(g.notpi, 2, mean)
apply(g.notpi, 2, sd)




## URM

nrows<-length(g.urm[,1])
ncols<-6 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==1) {
      figdata[i+5*nrows,1] <- g.urm$S8.alpha[i]
      figdata[i+5*nrows,2] <- j}
    if (j==2) {
      figdata[i+4*nrows,1] <- g.urm$S3.alpha[i]
      figdata[i+4*nrows,2] <- j}
    if (j==3) {
      figdata[i+3*nrows,1] <- g.urm$S2.alpha[i]
      figdata[i+3*nrows,2] <- j}
    if (j==4) {
      figdata[i+2*nrows,1] <- g.urm$S5.alpha[i]
      figdata[i+2*nrows,2] <- j}
    if (j==5) {
      figdata[i+1*nrows,1] <- g.urm$S6.alpha[i]
      figdata[i+1*nrows,2] <- j}
    if (j==6) {
      figdata[i,1] <- g.urm$S1.alpha[i]
      figdata[i,2] <- j}
  }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3", "4", "5", "6"),
  labels = c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy",
             "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates - URM Only',
    subtitle = 'Cohen\'s D',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())



## Non-Male


nrows<-length(g.notmale[,1])
ncols<-6 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==1) {
      figdata[i+5*nrows,1] <- g.notmale$S8.alpha[i]
      figdata[i+5*nrows,2] <- j}
    if (j==2) {
      figdata[i+4*nrows,1] <- g.notmale$S3.alpha[i]
      figdata[i+4*nrows,2] <- j}
    if (j==3) {
      figdata[i+3*nrows,1] <- g.notmale$S2.alpha[i]
      figdata[i+3*nrows,2] <- j}
    if (j==4) {
      figdata[i+2*nrows,1] <- g.notmale$S5.alpha[i]
      figdata[i+2*nrows,2] <- j}
    if (j==5) {
      figdata[i+1*nrows,1] <- g.notmale$S6.alpha[i]
      figdata[i+1*nrows,2] <- j}
    if (j==6) {
      figdata[i,1] <- g.notmale$S1.alpha[i]
      figdata[i,2] <- j}
  }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3", "4", "5", "6"),
  labels = c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy",
             "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates - Non-Male Only',
    subtitle = 'Cohen\'s D',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())


## Non-PI




nrows<-length(g.notpi[,1])
ncols<-6 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==1) {
      figdata[i+5*nrows,1] <- g.notpi$S8.alpha[i]
      figdata[i+5*nrows,2] <- j}
    if (j==2) {
      figdata[i+4*nrows,1] <- g.notpi$S3.alpha[i]
      figdata[i+4*nrows,2] <- j}
    if (j==3) {
      figdata[i+3*nrows,1] <- g.notpi$S2.alpha[i]
      figdata[i+3*nrows,2] <- j}
    if (j==4) {
      figdata[i+2*nrows,1] <- g.notpi$S5.alpha[i]
      figdata[i+2*nrows,2] <- j}
    if (j==5) {
      figdata[i+1*nrows,1] <- g.notpi$S6.alpha[i]
      figdata[i+1*nrows,2] <- j}
    if (j==6) {
      figdata[i,1] <- g.notpi$S1.alpha[i]
      figdata[i,2] <- j}
  }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3", "4", "5", "6"),
  labels = c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy",
             "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates - Not-PI Only',
    subtitle = 'Cohen\'s D',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())


## WM PI



nrows<-length(g.wmpi[,1])
ncols<-6 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==1) {
      figdata[i+5*nrows,1] <- g.wmpi$S8.alpha[i]
      figdata[i+5*nrows,2] <- j}
    if (j==2) {
      figdata[i+4*nrows,1] <- g.wmpi$S3.alpha[i]
      figdata[i+4*nrows,2] <- j}
    if (j==3) {
      figdata[i+3*nrows,1] <- g.wmpi$S2.alpha[i]
      figdata[i+3*nrows,2] <- j}
    if (j==4) {
      figdata[i+2*nrows,1] <- g.wmpi$S5.alpha[i]
      figdata[i+2*nrows,2] <- j}
    if (j==5) {
      figdata[i+1*nrows,1] <- g.wmpi$S6.alpha[i]
      figdata[i+1*nrows,2] <- j}
    if (j==6) {
      figdata[i,1] <- g.wmpi$S1.alpha[i]
      figdata[i,2] <- j}
  }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3", "4", "5", "6"),
  labels = c("Preserve Replication Materials", "Reasons for Data Management Policy", "Reasons for Authorship Policy",
             "Lab Disagreement", "Respectful Discussion", "Relevance of Ethics Discourse")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates - White, Male PI Only',
    subtitle = 'Cohen\'s D',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())






